{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "entitled-melbourne",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solving nonlinear systems of equations\n",
    "\n",
    "# To solve a system of nonlinear equations, we will use 'fsolve()' which requires\n",
    "# the `scipy.optimize' module.\n",
    "import numpy as np\n",
    "from scipy.optimize import fsolve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "comparable-light",
   "metadata": {},
   "outputs": [],
   "source": [
    "# First, we will solve a single nonlinear equations.  We start by defining\n",
    "# our equation as a function.\n",
    "def f(x):\n",
    "    return 3*x**3 - 2*x**2 + x - 7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "limiting-beaver",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-5"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Once the function is defined, we can evaluate it for different values of x\n",
    "# by:\n",
    "f(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "offensive-extreme",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can now use fsolve().  The first arguement is the function (i.e. the equation\n",
    "# that we want to solve, assumed to be equal to zero) and the second argument\n",
    "# is an initial guess at the solution.\n",
    "x = fsolve(f, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "focal-calcium",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.49175209])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The output of fsolve() is an array.\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "comfortable-command",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.4917520892393685"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# We can index the array to get the final solution.\n",
    "x[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "alone-tender",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1.2434497875801753e-14"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Let's check that the solution works.\n",
    "f(x[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "tender-memorabilia",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here's another example: Find y such that tan(e^(-2y)) = 1/y.\n",
    "def f(y):\n",
    "    return np.tan(np.exp(-2*y)) - 1/y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "alternate-powder",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.31441845309470834"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Here's the solution to f(y) = 0.\n",
    "y = fsolve(f, -0.1)\n",
    "y[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "silver-medline",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-2.220446049250313e-15"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Test the solution...\n",
    "f(y[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "enabling-values",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "y = -0.31441845309470834\n",
      "f(y) = -2.220446049250313e-15\n"
     ]
    }
   ],
   "source": [
    "# Note that there is another way to use fsolve such that it outputs just the \n",
    "# numerical solution (rather than inserting it into an array).\n",
    "y, = fsolve(f, -0.1)\n",
    "print('y =', y)\n",
    "print('f(y) =', f(y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "located-victory",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can also solve a system of nonlinear equations.  Below we define a system\n",
    "# of three equations with three unknown variables.\n",
    "def equations(vars):\n",
    "    x, y, z = vars\n",
    "    eq1 = 3*x*y + y - z -10\n",
    "    eq2 = x + x**2*y + z -12\n",
    "    eq3 = x -y -z + 2\n",
    "    return [eq1, eq2, eq3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "sustained-preference",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[-7, -14, 6]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# To evaluate the three equations at particular values of x, y, and z, we can\n",
    "# enter:\n",
    "equations([1, 0, -3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "coastal-reservoir",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here is the fsolve() statement with an initial guess at the solution\n",
    "x, y, z =  fsolve(equations, (1, 1, 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "bigger-screw",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x = 2.100953871762653 y = 1.6983245687167783 z = 2.4026293030458747\n"
     ]
    }
   ],
   "source": [
    "# Using this 'comma' notation, the get the three solutions directly.\n",
    "print('x =', x, 'y =', y, 'z =', z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "thorough-bottom",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "eq1 = 1.36e-10 eq2 = 2.76e-10 eq3 = 0\n"
     ]
    }
   ],
   "source": [
    "# Here, we check that all three of the equations actually evaluate to something\n",
    "# close to zero.  Note that, the 'SymPy' module is needed to use the 'N()' function. \n",
    "from sympy import *\n",
    "print('eq1 =', N(equations([x, y, z])[0], 3), 'eq2 =', N(equations([x, y, z])[1], 3),\\\n",
    "      'eq3 =', N(equations([x, y, z])[2], 3))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}